{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n",
      "70\n",
      "71\n",
      "72\n",
      "73\n",
      "74\n",
      "75\n",
      "76\n",
      "77\n",
      "78\n",
      "79\n",
      "80\n",
      "81\n",
      "82\n",
      "83\n",
      "84\n",
      "85\n",
      "86\n",
      "87\n",
      "88\n",
      "89\n",
      "90\n",
      "91\n",
      "92\n",
      "93\n",
      "94\n",
      "95\n",
      "96\n",
      "97\n",
      "98\n",
      "99\n",
      "100\n",
      "101\n",
      "102\n",
      "103\n",
      "104\n",
      "105\n",
      "106\n",
      "107\n",
      "108\n",
      "109\n",
      "110\n",
      "111\n",
      "112\n",
      "113\n",
      "114\n",
      "115\n",
      "116\n",
      "117\n",
      "118\n",
      "119\n",
      "120\n",
      "121\n",
      "122\n",
      "123\n",
      "124\n",
      "125\n",
      "126\n",
      "127\n",
      "128\n",
      "129\n",
      "130\n",
      "131\n",
      "132\n",
      "133\n",
      "134\n",
      "135\n",
      "136\n",
      "137\n",
      "138\n",
      "139\n",
      "140\n",
      "141\n",
      "142\n",
      "143\n",
      "144\n",
      "145\n",
      "146\n",
      "147\n",
      "148\n",
      "149\n",
      "150\n",
      "151\n",
      "152\n",
      "153\n",
      "154\n",
      "155\n",
      "156\n",
      "157\n",
      "158\n",
      "       Beijing_S Beijing_E  Beijing_I   Beijing_R Beijing_D   Beijing_H  \\\n",
      "0    2.15463e+07     0.059  0.0865333           0         0           0   \n",
      "1    2.15506e+07  0.131757   0.176054  0.00454243         0  0.00378205   \n",
      "2    2.15548e+07  0.218746   0.269829   0.0137436         0    0.012228   \n",
      "3    2.15591e+07  0.320566   0.369125   0.0277651         0   0.0262502   \n",
      "4    2.15634e+07  0.437945   0.475219   0.0468141         0   0.0467993   \n",
      "..           ...       ...        ...         ...       ...         ...   \n",
      "154  2.17244e+07    10.354    46.8995     233.145         0     261.532   \n",
      "155  2.17244e+07   10.2246    46.3323     234.376         0     262.196   \n",
      "156  2.17244e+07    10.097    45.7716     235.592         0     262.851   \n",
      "157  2.17244e+07   9.97123    45.2175     236.794         0     263.499   \n",
      "158  2.17244e+07   9.84714    44.6699      237.98         0     264.138   \n",
      "\n",
      "         Wuhan_S  Wuhan_E  Wuhan_I  Wuhan_R Wuhan_D  Wuhan_H  \n",
      "0    1.08886e+07  169.925  223.755  5.77428       1  50.6154  \n",
      "1    1.08843e+07  190.458  228.686  11.6448       1   61.508  \n",
      "2      1.088e+07  211.705  234.803  17.6425       1  73.7169  \n",
      "3    1.08756e+07   233.77  242.119  23.7984       1  87.2877  \n",
      "4    1.08713e+07  256.762  250.654  30.1438       1  102.273  \n",
      "..           ...      ...      ...      ...     ...      ...  \n",
      "154  1.06823e+07  692.458  4248.98    23254       1  27342.5  \n",
      "155  1.06823e+07  680.792  4181.84  23365.5       1  27386.9  \n",
      "156  1.06823e+07  669.351  4115.72  23475.3       1  27430.5  \n",
      "157  1.06823e+07   658.13  4050.61  23583.3       1  27473.4  \n",
      "158  1.06822e+07  647.123  3986.48  23689.6       1  27515.6  \n",
      "\n",
      "[159 rows x 12 columns]\n",
      "Running time: 2592.1657309532166 Seconds\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import time\n",
    "start=time.time()\n",
    "\n",
    "# 模拟天数 \n",
    "Tspan=53\n",
    "dt=1/3\n",
    "Delay=-3\n",
    "Tlock=14#封城时间 1月23日封城\n",
    "\n",
    "Tbeijinghome=15+Delay # 1月24日一级响应\n",
    "Tbeijingrelease=32+Delay  #逐渐恢复市内交通\n",
    "Twuhanhome=14+Delay # 武汉home\n",
    "Thubeirelease=80+Delay#湖北三月底恢复交通\n",
    "\n",
    "Trelease=140+Delay#湖北和北京交通三月初开通\n",
    "#城市属性\n",
    "Citys=['Beijing','Wuhan']\n",
    "Wuhanqu=['江岸区','江汉区','硚口区','汉阳区','武昌区','青山区',\n",
    "         '洪山区','蔡甸区','江夏区','黄陂区','新洲区','东西湖区','汉南区']\n",
    "Beijingqu=['密云区','延庆区','朝阳区','丰台区','石景山区','海淀区','门头沟区', \n",
    "           '房山区','通州区','顺义区','昌平区','大兴区','怀柔区','平谷区','东城区','西城区']\n",
    "Wuhan_renkou=10000*np.array([96.24,72.96,86.85,65.27,127.63,52.88,163.75,\n",
    "                       73.50,91.37,98.83,90.21,56.25,13.55])\n",
    "\n",
    "WUHAN=pd.DataFrame(data=Wuhan_renkou,index=Wuhanqu,columns=['Renkou'])\n",
    "Beijing_renkou=10000*np.array([49.5,34.8,360.5,210.5,59.0,335.8,33.1,118.8,157.8,116.9,210.8,179.6,\n",
    "                        41.4,45.6,82.2,117.9])\n",
    "BEIJING=pd.DataFrame(data=Beijing_renkou,index=Beijingqu,columns=['Renkou'])\n",
    "Patches={'Beijing':Beijingqu,'Wuhan':Wuhanqu}\n",
    "\n",
    "#Wuhancases=np.array([258,363,425,495,572,618,698,1590,1905,2261,2639,3215,4109,5142,6384,8351,10117,11618,13603,14982,16902,18454,19558,...\n",
    "   #20630+12364,21960+14031,37914,39462,41152,42752,44412,45027,45346,45660])\n",
    "\n",
    "Atr=pd.DataFrame(index=Citys)#城市属性\n",
    "Atr['Pnum']=np.zeros(len(Citys))\n",
    "for i in Citys:#添加城市属性Patch num区的数目\n",
    "    Atr.loc[i,'Pnum']=len(Patches[i])\n",
    "Atr['Province']=['Beijing','Hubei'] #添加城市省份\n",
    "\n",
    "\n",
    "#c*w*h个index\n",
    "INDEX=[]\n",
    "for c in Citys:\n",
    "    for w in Patches[c]:\n",
    "        for h in Patches[c]:\n",
    "            INDEX.append(c+h+w)\n",
    "            \n",
    "            \n",
    "# 输入的每个index的参数属性\n",
    "PARA=pd.DataFrame(index=INDEX,columns=['Cityname','S_initial','E_initial',\n",
    "                           'I_initial','R_initial','D_initial','H_initial'])\n",
    "for c in Citys:#更新城市名字\n",
    "    for i in PARA.index:\n",
    "        if c in i:\n",
    "            PARA.loc[i,'Cityname']=c\n",
    "PARA['S_initial']=np.ones(len(INDEX))\n",
    "c='Wuhan'\n",
    "for h in Patches[c]:\n",
    "    for w in Patches[c]:\n",
    "        if h==w:\n",
    "            PARA.loc[c+h+w,'S_initial']=WUHAN.loc[h,'Renkou']/3+WUHAN.loc[h,'Renkou']**2/sum(Wuhan_renkou)*2/3\n",
    "        else:\n",
    "            PARA.loc[c+h+w,'S_initial']=2/3*WUHAN.loc[h,'Renkou']*WUHAN.loc[w,'Renkou']/sum(Wuhan_renkou)\n",
    "c='Beijing'\n",
    "for h in Patches[c]:\n",
    "    for w in Patches[c]:\n",
    "        if h==w:\n",
    "            PARA.loc[c+h+w,'S_initial']=BEIJING.loc[h,'Renkou']/3+BEIJING.loc[h,'Renkou']**2/sum(Beijing_renkou)*2/3\n",
    "        else:\n",
    "            PARA.loc[c+h+w,'S_initial']=2/3*BEIJING.loc[h,'Renkou']*BEIJING.loc[w,'Renkou']/sum(Beijing_renkou)    \n",
    "'''AAA=0\n",
    "c='Wuhan'\n",
    "for h in Patches[c]:\n",
    "    for w in Patches[c]:\n",
    "        AAA=AAA+PARA.loc[c+h+w,'S_initial']\n",
    "print(AAA)'''\n",
    "PARA['E_initial']=np.zeros(len(INDEX))\n",
    "PARA.loc['Wuhan江汉区江汉区','E_initial']=150\n",
    "PARA['I_initial']=np.zeros(len(INDEX))\n",
    "PARA.loc['Wuhan江汉区江汉区','I_initial']=220\n",
    "PARA['R_initial']=np.zeros(len(INDEX))\n",
    "PARA['D_initial']=np.zeros(len(INDEX))\n",
    "PARA.loc['Wuhan江汉区江汉区','D_initial']=1\n",
    "PARA['H_initial']=np.zeros(len(INDEX))\n",
    "PARA.loc['Wuhan江汉区江汉区','H_initial']=41\n",
    "#print(PARA)\n",
    "\n",
    "\n",
    "# 输出\n",
    "#各个地方的S^{c,w,h}\n",
    "COLUMNS=[] #每列\n",
    "for i in INDEX:\n",
    "    A=[i+'_S',i+'_E',i+'_I',i+'_R',i+'_D',i+'_H']\n",
    "    COLUMNS.extend(A)\n",
    "#print(COLUMNS)\n",
    "Solution=pd.DataFrame(columns=COLUMNS)\n",
    "\n",
    "\n",
    "#print(Solution)\n",
    "def Control(c,i,t): #city, location, time, initial time is Jan 9th,24:00 OK\n",
    "    cnh=20+Delay#控制措施产生效果\n",
    "    cnhk=0.14#效果\n",
    "    cnhday=5#花了多少天\n",
    "    cwh=20+Delay#控制措施产生效果\n",
    "    cwhday=8#花了多少天\n",
    "    cwhk=0.083#效果\n",
    "    def cNH(t):#湖北以外控制\n",
    "        if t<=cnh: #1月30 号之后\n",
    "            c=1\n",
    "        elif t>cnh+cnhday:\n",
    "            c=cnhk  #可调\n",
    "        else:\n",
    "            c=np.exp(np.log(cnhk)*(t-cnh)/cnhday)\n",
    "        return c\n",
    "    def cH(t):#湖北控制\n",
    "        if t<=cwh: #2月3号之后\n",
    "            c=1\n",
    "        elif t>cwh+cwhday:\n",
    "            c=cwhk  #可调\n",
    "        else:\n",
    "            c=np.exp(np.log(cwhk)*(t-cwh)/cwhday)\n",
    "        return c\n",
    "    if Atr.loc[c,'Province']=='Hubei':\n",
    "        control=cH(t)\n",
    "    else:\n",
    "        control=cNH(t)\n",
    "    return control\n",
    "#print(Control('Wuhan',16,26))\n",
    "def Betae(c,i,t):#Exposed 感染率， i表示location, t时间 OK\n",
    "    if Atr.loc[c,'Province']=='Hubei':\n",
    "        betae=0.24\n",
    "    else:\n",
    "        betae=0.20\n",
    "    return betae\n",
    "#print(Betae('Beijing',10))\n",
    "def Betai(c,i,t):#Infected 感染率， i表示location, t时间 OK\n",
    "    if Atr.loc[c,'Province']=='Hubei':\n",
    "        betai=0.24\n",
    "    else:\n",
    "        betai=0.20\n",
    "    return betai\n",
    "def Prob(c,w,h,i,t): # c city, w work place, h, home place, i, location, t time OK\n",
    "    if c=='Beijing':\n",
    "        if t<=Tbeijinghome:   \n",
    "            prob=1/Atr.loc[c,'Pnum']\n",
    "        elif t>Tbeijinghome and t<Tbeijingrelease:\n",
    "            if i==h:\n",
    "                prob=1\n",
    "            else:\n",
    "                prob=0\n",
    "        else:\n",
    "            if (np.floor(t)-Tbeijingrelease+1)%7==5 or (np.floor(t)-Tbeijingrelease+1)%7==6:#双休日\n",
    "                if t-np.floor(t)>0 and t-np.floor(t)<=1/3+0.001:#白天\n",
    "                    prob=1/Atr.loc[c,'Pnum']\n",
    "                else:\n",
    "                    if i==h:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "            else:\n",
    "                if t-np.floor(t)>0 and t-np.floor(t)<=1/3+0.001:#白天\n",
    "                    if i==w:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "                else:\n",
    "                    if i==h:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "    if c=='Wuhan':\n",
    "        if t<=Twuhanhome:   \n",
    "            prob=1/Atr.loc[c,'Pnum']\n",
    "        elif t>Twuhanhome and t<Thubeirelease:\n",
    "            if i==h:\n",
    "                prob=1\n",
    "            else:\n",
    "                prob=0\n",
    "        else:\n",
    "            if (np.floor(t)-Thubeirelease+1)%7==5 or (np.floor(t)-Thubeirelease+1)%7==6:#双休日\n",
    "                if t-np.floor(t)>0 and t-np.floor(t)<=1/3+0.001:#白天\n",
    "                    prob=1/Atr.loc[c,'Pnum']\n",
    "                else:\n",
    "                    if i==h:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "            else:#工作日\n",
    "                if t-np.floor(t)>0 and t-np.floor(t)<=1/3+0.001:#白天\n",
    "                    if i==w:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "                else:\n",
    "                    if i==h:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0            \n",
    "    return prob\n",
    "'''for i in range(180):\n",
    "    print(Prob('Wuhan','汉南区','江汉区','汉南区',i/3))'''\n",
    "def Diffusion(a,b,t): # a city to b city, distance induced diffusion rate OK\n",
    "    if t>Tlock and t<Trelease:\n",
    "        diffusion=0\n",
    "    else:\n",
    "        diffusion=0.00118\n",
    "    return diffusion\n",
    "def Migration(a,b,t):#a city to b city,migration probability OK\n",
    "    if a=='Wuhan' and b=='Beijing':\n",
    "        migration=1\n",
    "    else:\n",
    "        migration=0\n",
    "    return migration\n",
    "print(Migration('Wuhan','Beijing',2))\n",
    "def Qprob(c,w,h,t):# c city, w work place, h, home place, t time OK\n",
    "    qprob=1/Atr.loc[c,'Pnum']**2\n",
    "    return qprob\n",
    "#print(Qprob('Beijing',1,1,1))\n",
    "def Gamma(c,i,t): # c city, i patch, t time OK\n",
    "    y=1/12.7\n",
    "    return y\n",
    "def Alpha(c,i,t): # c city, i patch, t time OK\n",
    "    a=0;\n",
    "    return a\n",
    "def Sigma(): #OK\n",
    "    s=1/5.2\n",
    "    return s\n",
    "\n",
    "#定义中间变量\n",
    "COLUMNS1=[] #每列\n",
    "for i in INDEX:\n",
    "    A=[i+'_S1',i+'_E1',i+'_I1',i+'_R1',i+'_D1',i+'_H1']\n",
    "    COLUMNS1.extend(A)\n",
    "COLUMNS2=[] #每列\n",
    "for i in INDEX:\n",
    "    A=[i+'_S2',i+'_E2',i+'_I2',i+'_R2',i+'_D2',i+'_H2']\n",
    "    COLUMNS2.extend(A)\n",
    "#print(COLUMNS1)\n",
    "Midpar1=pd.DataFrame(data=np.zeros((1,len(COLUMNS1))),columns=COLUMNS1)\n",
    "Midpar2=pd.DataFrame(data=np.zeros((1,len(COLUMNS2))),columns=COLUMNS2)\n",
    "#print(type(Midpar2.iloc[:,0]))\n",
    "#代入初值\n",
    "for i in INDEX:\n",
    "    for j in ['S','E','I','R','D','H']:\n",
    "        Midpar1[i+'_'+j+'1']=PARA.loc[i,j+'_initial']\n",
    "#print(Midpar1)\n",
    "for i in range(int(Tspan/dt)):\n",
    "    print(i)\n",
    "    for c in Citys:\n",
    "        #计算S_{i}^{c} E_{i}^{c} I R_{i}^{c}\n",
    "        Nlocation=pd.DataFrame(data=np.zeros((len(Patches[c]),5)),\n",
    "                          columns=['Ssum','Esum','Isum','Rsum','Nsum'],index=Patches[c])\n",
    "        #print(type(Nlocation.loc[patch,'Ssum']))\n",
    "        for patch in Patches[c]:\n",
    "            for hh in Patches[c]:\n",
    "                for ww in Patches[c]:\n",
    "                    Suma=Nlocation.loc[patch,'Ssum']+Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_S1']\n",
    "                    #print(Suma)\n",
    "                    #print(type(Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_S1']))\n",
    "                    Sumb=Nlocation.loc[patch,'Esum']+Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_E1']\n",
    "                    Sumc=Nlocation.loc[patch,'Isum']+Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_I1']\n",
    "                    Sumd=Nlocation.loc[patch,'Rsum']+Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_R1']\n",
    "                    Nlocation.loc[patch,'Ssum']=Suma\n",
    "                    Nlocation.loc[patch,'Esum']=Sumb\n",
    "                    Nlocation.loc[patch,'Isum']=Sumc\n",
    "                    Nlocation.loc[patch,'Rsum']=Sumd\n",
    "                    Nlocation.loc[patch,'Nsum']=Suma+Sumb+Sumc+Sumd\n",
    "        Ntotal=sum(Nlocation['Nsum'].values)\n",
    "        #print(Ntotal) OK\n",
    "        #计算扩散项:迁出 和迁入 \n",
    "        imigration_rate=0\n",
    "        emigration_num_S=0\n",
    "        emigration_num_E=0\n",
    "        emigration_num_I=0\n",
    "        emigration_num_R=0\n",
    "        for b in Citys:\n",
    "            imigration_rate=imigration_rate+Migration(c,b,i*dt)*Diffusion(c,b,i*dt)\n",
    "            #计算 S^{b}\n",
    "            SSUM=0\n",
    "            ESUM=0\n",
    "            ISUM=0\n",
    "            RSUM=0\n",
    "            for hb in Patches[b]:\n",
    "                for wb in Patches[b]:\n",
    "                    SSUM=SSUM+Midpar1.loc[0,b+hb+wb+'_S1']\n",
    "                    #print(SSUM)\n",
    "                    ESUM=ESUM+Midpar1.loc[0,b+hb+wb+'_E1']\n",
    "                    ISUM=ISUM+Midpar1.loc[0,b+hb+wb+'_I1']\n",
    "                    RSUM=RSUM+Midpar1.loc[0,b+hb+wb+'_R1']\n",
    "            emigration_num_S=emigration_num_S+Migration(b,c,i*dt)*Diffusion(b,c,i*dt)*SSUM\n",
    "            emigration_num_E=emigration_num_E+Migration(b,c,i*dt)*Diffusion(b,c,i*dt)*ESUM\n",
    "            emigration_num_I=emigration_num_I+Migration(b,c,i*dt)*Diffusion(b,c,i*dt)*ISUM\n",
    "            emigration_num_R=emigration_num_R+Migration(b,c,i*dt)*Diffusion(b,c,i*dt)*RSUM\n",
    "        for h in Patches[c]:\n",
    "            for w in Patches[c]:       \n",
    "                Infection=0\n",
    "                Recover=0\n",
    "                Dead=0\n",
    "                for patch in Patches[c]:\n",
    "                    I_sum=Nlocation.loc[patch,'Isum']\n",
    "                    E_sum=Nlocation.loc[patch,'Esum']\n",
    "                    N_sum=Nlocation.loc[patch,'Nsum']\n",
    "                    Infection=Infection+Prob(c,w,h,patch,i*dt)*Control(c,patch,i*dt)*(Betai(c,patch,i*dt)*I_sum/N_sum+Betae(c,patch,i*dt)*E_sum/N_sum)\n",
    "                    Recover=Recover+Prob(c,w,h,patch,i*dt)*Gamma(c,patch,i*dt)\n",
    "                    Dead=Dead+Prob(c,w,h,patch,i*dt)*Alpha(c,patch,i*dt)\n",
    "                #print(Recover) OK\n",
    "                qprob=(Midpar1.loc[0,c+h+w+'_S1']+Midpar1.loc[0,c+h+w+'_E1']+Midpar1.loc[0,c+h+w+'_I1']+Midpar1.loc[0,c+h+w+'_R1'])/Ntotal\n",
    "                Midpar2.loc[0,c+h+w+'_S2']=Midpar1.loc[0,c+h+w+'_S1']+dt*(-Midpar1.loc[0,c+h+w+'_S1']*Infection-imigration_rate*Midpar1.loc[0,c+h+w+'_S1']+qprob*emigration_num_S)\n",
    "                Midpar2.loc[0,c+h+w+'_E2']=Midpar1.loc[0,c+h+w+'_E1']+dt*(Midpar1.loc[0,c+h+w+'_S1']*Infection-Sigma()*Midpar1.loc[0,c+h+w+'_E1']-imigration_rate*Midpar1.loc[0,c+h+w+'_E1']+qprob*emigration_num_E)\n",
    "                Midpar2.loc[0,c+h+w+'_I2']=Midpar1.loc[0,c+h+w+'_I1']+dt*(Sigma()*Midpar1.loc[0,c+h+w+'_E1']-(Recover+Dead)*Midpar1.loc[0,c+h+w+'_I1']-imigration_rate*Midpar1.loc[0,c+h+w+'_I1']+qprob*emigration_num_I)\n",
    "                Midpar2.loc[0,c+h+w+'_R2']=Midpar1.loc[0,c+h+w+'_R1']+dt*(Recover*Midpar1.loc[0,c+h+w+'_I1']-imigration_rate*Midpar1.loc[0,c+h+w+'_R1']+qprob*emigration_num_R)\n",
    "                Midpar2.loc[0,c+h+w+'_D2']=Midpar1.loc[0,c+h+w+'_D1']+dt*Dead*Midpar1.loc[0,c+h+w+'_I1']\n",
    "                Midpar2.loc[0,c+h+w+'_H2']=Midpar1.loc[0,c+h+w+'_H1']+dt*Sigma()*Midpar1.loc[0,c+h+w+'_E1']\n",
    "                Midpar1.loc[0,c+h+w+'_S1']=Midpar2.loc[0,c+h+w+'_S2']\n",
    "                Midpar1.loc[0,c+h+w+'_E1']=Midpar2.loc[0,c+h+w+'_E2']\n",
    "                Midpar1.loc[0,c+h+w+'_I1']=Midpar2.loc[0,c+h+w+'_I2']\n",
    "                Midpar1.loc[0,c+h+w+'_R1']=Midpar2.loc[0,c+h+w+'_R2']\n",
    "                Midpar1.loc[0,c+h+w+'_D1']=Midpar2.loc[0,c+h+w+'_D2']\n",
    "                Midpar1.loc[0,c+h+w+'_H1']=Midpar2.loc[0,c+h+w+'_H2']\n",
    "                Solution.loc[i,c+h+w+'_S']=Midpar2.loc[0,c+h+w+'_S2']\n",
    "                Solution.loc[i,c+h+w+'_E']=Midpar2.loc[0,c+h+w+'_E2']\n",
    "                Solution.loc[i,c+h+w+'_I']=Midpar2.loc[0,c+h+w+'_I2']\n",
    "                Solution.loc[i,c+h+w+'_R']=Midpar2.loc[0,c+h+w+'_R2']\n",
    "                Solution.loc[i,c+h+w+'_D']=Midpar2.loc[0,c+h+w+'_D2']\n",
    "                Solution.loc[i,c+h+w+'_H']=Midpar2.loc[0,c+h+w+'_H2']\n",
    "#S^{c,h}\n",
    "#各区感染人口\n",
    "INDEX_Distinct=[]\n",
    "for  c in Citys:\n",
    "    for h in Patches[c]:\n",
    "        INDEX_Distinct.append(c+h)\n",
    "COLUMNS_Distinct=[] #每列\n",
    "for i in INDEX_Distinct:\n",
    "    A=[i+'_S',i+'_E',i+'_I',i+'_R',i+'_D',i+'_H']\n",
    "    COLUMNS_Distinct.extend(A)\n",
    "Solution_Distinct=pd.DataFrame(columns=COLUMNS_Distinct)\n",
    "for i in range(int(Tspan/dt)):\n",
    "    for c in Citys:\n",
    "        for h in Patches[c]:\n",
    "            S_ch=0\n",
    "            E_ch=0\n",
    "            I_ch=0\n",
    "            R_ch=0\n",
    "            D_ch=0\n",
    "            H_ch=0\n",
    "            for w in Patches[c]:\n",
    "                S_ch=S_ch+Solution.loc[i,c+h+w+'_S']\n",
    "                E_ch=E_ch+Solution.loc[i,c+h+w+'_E']\n",
    "                I_ch=I_ch+Solution.loc[i,c+h+w+'_I']\n",
    "                R_ch=R_ch+Solution.loc[i,c+h+w+'_R']\n",
    "                D_ch=D_ch+Solution.loc[i,c+h+w+'_D']\n",
    "                H_ch=H_ch+Solution.loc[i,c+h+w+'_H']\n",
    "            Solution_Distinct.loc[i,c+h+'_S']=S_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_E']=E_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_I']=I_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_R']=R_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_D']=D_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_H']=H_ch\n",
    "#各市感染人口\n",
    "\n",
    "COLUMNS_City=[] #每列\n",
    "for i in Citys:\n",
    "    A=[i+'_S',i+'_E',i+'_I',i+'_R',i+'_D',i+'_H']\n",
    "    COLUMNS_City.extend(A)\n",
    "Solution_City=pd.DataFrame(columns=COLUMNS_City)           \n",
    "for i in range(int(Tspan/dt)):\n",
    "    for c in Citys:\n",
    "        S_c=0\n",
    "        E_c=0\n",
    "        I_c=0\n",
    "        R_c=0\n",
    "        D_c=0\n",
    "        H_c=0        \n",
    "        for h in Patches[c]:\n",
    "            S_c=S_c+Solution_Distinct.loc[i,c+h+'_S']\n",
    "            E_c=E_c+Solution_Distinct.loc[i,c+h+'_E']\n",
    "            I_c=I_c+Solution_Distinct.loc[i,c+h+'_I']\n",
    "            R_c=R_c+Solution_Distinct.loc[i,c+h+'_R']\n",
    "            D_c=D_c+Solution_Distinct.loc[i,c+h+'_D']\n",
    "            H_c=H_c+Solution_Distinct.loc[i,c+h+'_H']\n",
    "        Solution_City.loc[i,c+'_S']=S_c\n",
    "        Solution_City.loc[i,c+'_E']=E_c\n",
    "        Solution_City.loc[i,c+'_I']=I_c\n",
    "        Solution_City.loc[i,c+'_R']=R_c\n",
    "        Solution_City.loc[i,c+'_D']=D_c\n",
    "        Solution_City.loc[i,c+'_H']=H_c\n",
    "            \n",
    "print(Solution_City)\n",
    "end=time.time()\n",
    "print('Running time: %s Seconds'%(end-start))\n",
    "\n",
    "Solution_City.to_csv('CityBefore3City202003302.csv',index=0)\n",
    "Solution_Distinct.to_csv('CityBefore3Distinct202003302.csv',index=0)\n",
    "Solution.to_csv('CityBefore3Detail202003302.csv',index=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
